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ABSTRACT 

We investigate the acceleration of light particles in perpendicular shocks for plasmas 
consisting of a mixture of leptonic and hadronic particles. Starting from the full set of 
conservation equations for the mixed plasma constituents, we generalize the magnet o- 
hydrodynamical jump conditions for a multi-component plasma, including information 
about the specific adiabatic constants for the different species. The impact of devia- 
tions from the standard model of an ideal gas is compared in theory and particle-in-cell 
simulations, showing that the standard-MHD model is a good approximation. The simu- 
lations of shocks in electron-positron-ion plasmas are for the first time multi-dimensional, 
transverse effects are small in this configuration and ID simulations are a good repre- 
sentation if the initial magnetization is chosen high. ID runs with a mass ratio of 1836 
are performed, which identify the Larmor frequency uj c i as the dominant frequency that 
determines the shock physics in mixed component plasmas. The maximum energy in the 
non-thermal tail of the particle spectra evolves in time according to a power-law oc t a 
with a in the range 1/3 < a < 1, depending on the initial parameters. A connection 



is made with transport theoretical models by Drury (1983) and Gargate & Spitkovsky 



(2011 ), which predict an acceleration time oc 7 and the theory for small wavelength scat- 
tering by Kirk & Reville (2010), which predicts a behavior rather as oc 7 2 . Furthermore, 



we compare different magnetic field orientations with Bo inside and out of the plane, 
observing qualitatively different particle spectra than in pure electron-ion shocks. 

Subject headings: acceleration of particles, equation of state, ISM: kinematics and dy- 
namics, shock waves 



1. Introduction 



Shock acceleration has received considerable attention in recent years, due to the possibility 
of accelerating charged particles to very high energies. The investigation of shocks in pair plasmas 
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has been motivated by Arons (1983), arguing that these plasma constituents are dominant in some 



astrophysical scenarios. These scenarios are also convenient for numerical simulations (e. g. Chang et 
al. (2008), Spitkovsky (2008b), Nishikawa et al. ( 2009| )), due to the fact that numerical simulations 
including heavier particles are more demanding due to the disparity of typical length and time 
scales. Simulations of electron-ion shocks are mostly performed with a reduced mass ratio (e. g. 



Hoshino & Shimada (2002), Hededal et al. 


(2004 


), 


Spitkovsky 


(2008a 


), 



fc Takabe] ( |2010[ )). Recently, |Haugb0lle| ( |2010[ ) studied the full development and relaxation process 
of an electron-ion shock in a three-dimensional simulation, using a mass ratio rrii/m e = 16, in a 



two-dimensional spatial configuration, Sironi & Spitkovsky (2011) studied electron-ion shocks for 
the first time with mass ratios rrii/m e = 1000. 



Spitkovsky (2008a) showed that electron-ion shocks behave similarly as electron-positron shocks 
on large time scales, because the particle rest mass becomes negligible in comparison with the rela- 
tivistic mass once the stage of full downstream thermalization has been obtained. This is supported 
by the fact that both particle components show comparable energy spectra ( Martins et al.|20 09), and 
facilitates the comparison with theoretical models, as the standard jump conditions for a single-fluid, 



derived by Blandford & McKee (1976), can be applied. But this is true only for initially unmag- 
netized or quasi-parallel shocks. The investigation of strongly magnetized perpendicular shocks, 
that we perform in this paper, shows a different picture. The compression ratio is significantly 
increased in the presence of a heavy particle component and the shock front propagates at a lower 
velocity. The results are in good agreement with our analytical derivation of the jump conditions 
for perpendicular shocks in a multiple species plasma. For this analysis, the only assumption we 
make is that the downstream density profile is similar for all particle components, which happens 
after a few uj~^ even if the ions have not thermalized yet. 

There still exists a gap between analytical acceleration models and numerical simulations or 
observation data of perpendicular shocks, as very large amplitudes for the magnetic turbulence 
are needed to enable multiple scatterings of the particles in the shock and to form the observed 
energy spectra ( Arons||2007 ). Since the acceleration process cannot be explained by a simple model, 



Arons (2007) suggested a mixture of diffusive Fermi acceleration and an additional acceleration 
process, where heavy ions play a major role. The latter process provides a mechanism to produce a 
broad energy range in plasmas, where pairs dominate by number and where ions are energetically 
dominant, explaining the observed range of optical, X- and 7-radiation in Pulsar Wind Nebulae, 
but still the spectra are not all in agreement with observations. A characteristic of pure electron- 
positron perpendicular shocks is, that no evidence has been found for the existence of a non-thermal 



population (e.g. Langdon et al. (1988), Gallant et al. (1992)). This population appears only if the 



initial magnetic field has an oblique structure (Sironi & Spitkovsky 2009) or in the presence of 



an ion population (Hoshino et al. (1992), Amato & Arons ( 2006| )). In the latter case, Hoshino & 
Arons (1991) found that the light plasma species gains energy from the heavy species due to the 



synchrotron maser instability. The gyrating ions emit magnetosonic waves, which are absorbed 
preferentially by positrons, accelerating them to non-thermal energies. A ratio n^m^/n e m e > 10 
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is necessary to achieve efficient acceleration, as demonstrated by Amato & Arons (2006) in a ID 



simulation with rrii/m e — 100. The left-handed orientation of the waves, facilitates the energization 
of the positrons, which is why the electron spectrum in such a configuration was not observed to 



reach the same level as that of the positrons. Hoshino et al. (1992) suggested that a realistic mass 



ratio rrii/m e = 1836, will have the same accelerating effect on the electrons. We confirm that the 
electron tail is stronger in this case. However, we observe the acceleration efficiency to be not only 
a function of the ion mass but also the total magnetization. 

The temporal evolution of the maximum energy is investigated for different ion to electron 
density ratios and it is found to be consistent with the acceleration model due to multiple scattering 



in small wavelength turbulence (Kirk & Reville 2010). 



This paper is structured as follows. The physical scenario is described in Section [2] and the 
jump conditions are presented for a perpendicular shock in a plasma consisting of mixed particle 
constituents with different energy spectra. In Section [3] the simulation results are compared with 
the theory from the previous section. We discuss first the differences in the particle spectra and 
their effects on the jump conditions, where we vary the initial ion kinetic energy ratio for a constant 
magnetization, which leads to the same jump conditions in the standard MHD model. After, we 
discuss the advanced model for a wide range of parameters with decreasing magnetization. Finally, 
the effect of the magnetic field orientation is discussed briefly and the main results of the simulations 
are discussed and summarized in Section [U 



2. Theory for mixed particle components 

We investigate the interaction of two counterstreaming beams, where each stream is charge- 
neutral and consisting of a mixture of electrons, positrons and ions, in a constant perpendicular 
magnetic field. This leads to the formation of two shocks, each propagating in the opposite direction 
of the incoming upstream beam, away from the interaction region. The following analytical model 
describes the quasi-steady state after the shock is formed and our discussions throughout this paper 
are limited to the description of the shock propagating to the right-hand-side (see Figure [I]) . We 



adopt the syntax of Zhang & Kobayashi (2005) denoting quantities measured in their rest frame with 
a single index Qi and quantities measured in the rest frame j by Qij with z, j = 1, 2, s, denoting 
the upstream, downstream and shock frame, respectively. The additional index separated by a 
comma Qi, a , Qij, a with a = e— , e+, p+ specifies the species, electron, positron, ion, respectively. 
As the theory is compared with the simulation results in the following section, the calculations are 
performed in the simulation frame, which coincides with the downstream frame. In this frame, 
the different species can be treated equally with = ^2 = and 72 = 1. Moreover, all three 
components have the same upstream velocity /?i2 = vu/c < and Lorentz factor 712 and the shock 
propagates with /3 := /3 S 2 = Vs2 /c > and its associated Lorentz factor 7 := 7^2- 



Because of the charge-neutrality condition, the initial upstream densities of positrons and ions 
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Fig. 1. — Shock setup and definition of the right-hand-side (red box). 



sum up to the electron density nu := ni2,e- = ^i2,e+ + ^I2,p+- The simulations show that the 
downstream densities of the three species behave similarly after a few ou~J-, even if the ions have 
not thermalized yet, therefore we assume r^e-A^e- — ^2,e+/^i2,e+ + ^2,p+/^i2,p+- With these 
assumptions the jump conditions can be derived in a similar way as for a single fluid. The detailed 
derivation of the jump conditions is presented in the Appendix. For the remainder of the paper, we 
will use the shock speed, as a function of the key parameters, given by 
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and the approximation of Equation ([!]) for highly relativistic upstream Lorentz factors 712 ^> 1 
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Equation ^ reduces to Equation (16) of |Gallant et al. (1992) in the limit of equal downstream 
spectra, where in Equation ^ it is important to keep in mind that the effective magnetization g\ — 
a to ti defined in Equation ( |A7| ), is considered, containing the contributions of all particle components. 
The jump conditions are determined by the parameter 01, and the downstream adiabatic constant 
r e _. The dependence of the shock speed, defined by Equation ([2]), on the total magnetization and 
the initial ion to electron kinetic energy fraction m p ni2, p + / (m e ni2, e -) is demonstrated in Figure |2j 
For the sake of simplicity equal particle spectra are chosen with an adiabatic constant T e _ = 3/2. 
The shock speed is increased if the total magnetization is increased and decreased with increasing 
initial ion kinetic energy ratio. 



The jump conditions for the density and magnetic field ratios obtained from Equations (Al) 
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Fig. 2. — Shock speed f3 against (a) the total magnetization g\ and (b) the fraction of upstream 
kinetic energy carried by hadrons U p /U e — m p ni2, p + / (m e ni2, e -) with m p /m e = 100 and a e = 2 for 
T e _ = 1.5 (solid), 1.7 (large dashed), 1.9 (small dashed), assuming equal particle spectra. 



and (A2) are equal and given by 
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(3) 



2.1. The role of the adiabatic constant in the shock properties 



One of the objectives of this work is to determine the impact of the real particle distributions 
on the jump conditions. For this, the downstream adiabatic constants and pressure densities have 
to be determined from the simulation data. In the previous section, the adiabatic constant has been 
defined for each species as a relation between the energy, pressure and spatial densities, which are 
defined by 
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in a two-dimensional geometry. After the particle distribution has been determined from the sim- 
ulation data, the pressure densities and adiabatic constants T a — 1 + P2,a/(^2,a — ^2,a^aC 2 ) can be 
determined. The distribution functions are found to be fitted well by a Maxwellian for low energies 
plus a high-energy power-law tail and an exponential cut-off (e.g. Spitkovsky||2008b ) 

-\dn _ ^ r _ /A _] , n _ -(l+a) 
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exp [-7/A7] + C27' 



min{l, exp [-(7 - j cut ) / 'Aj cut ]} 
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with C2 — for 7 < jmin. An analytical expression of the densities Q-Q is provided in Stockem 
et al. (2011). A parameter study of Equation ^ is presented in Figure [3] showing the variation of 
the shock speed with the downstream adiabatic constant for a particular initial magnetization. If 
the magnetic field is strong, the impact of the change in the adiabatic constant, which is determined 
by the shape of the distribution function, is low, unlike what has been observed for unmagnetized 



scenarios m IStockem et all (12011}. In the unmagnetized case, the deviation in the shock speed will 
be 20%, whereas it is just 12% for <j\ = 0.2 or 5% for o\ = 1. 




Fig. 3. — Deviations in the shock speed A/3 = (3o — with /?o the shock speed for T e _ = 3/2, 
according to variations in the adiabatic constant Ar e _ = T e _ — 3/2 for magnetizations g\ = 
(solid), 0.2 (large dashed), 1 (small dashed), 10 (dotted). 



3. Simulations of highly magnetized perpendicular shocks 

To study the effect of mixed plasma components on the jump conditions and to test the theory 
developed in the previous section, we use ID and 2D particle-in-cell (PIC) simulations, which we 
perform with the kinetic PIC code OSIRIS ( |Fonseca et~aT1|2002[ |2008[ ). We have found that the 



setup described in Figure [T] is more appropriate for the numerical study of these shocks, since it 
avoids boundary condition issues. In this work, we employ a setup reproducing the model of Figure 
[TJ using two counterstreaming beams, so that two shocks are formed, propagating in opposite 
directions, away from the interaction region. The background magnetic field is constant in time and 
changes its sign according to the sign of the upstream velocities ±/3i2 of the opposite beams. The 
motional electric field E12 = x B12 is thus constant over the entire simulation box. Because 

of the symmetry of the formation of the two shocks, we limit our discussions to the right-hand side 
of the simulation box (see Figure [T]). 
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3.1. Varying the ion kinetic energy for a constant magnetization 

By adjusting the ion mass and density ratios, different ion kinetic energy ratios U p /U e — 
ni2, P +m p /ni2,e-me can lead to the same total magnetization o\ — a e /[2+(m p /m e — l)ni2,p+/ni2,e-] 3 
with the same jump conditions in the standard MHD model. In this section we discuss the differences 
in the spectra and their effects on the jump conditions. Due to the constraint < ni2,p+/^i2,e- < 1? 
the difference between the lowest and largest value of the ion kinetic energy ratio for constant u\ 
and a e is limited to A(U p /U e ) = 1. The total magnetization is chosen high enough (a\ > 10 -3 
(Spitkovsky 2005| )) to suppress the Weibel instability so that it can be excluded as the dominant 



driver for shock formation. A discussion of the role of the Weibel instability in baryon-loaded 



plasmas is presented in Fiore et al. (2006). 



We performed two sets of simulations for magnetizations u\ — 0.345 with ion kinetic energy 
ratios U p /U e — 4.0 and 4.5 and u\ — 0.145 with U p /U e — 12.0 and 12.5. The details are given 
in Fig. [ij Global parameters are the time step At = 0.2214 a;^ 1 = 3.9 ps x \/ n i2,e-[ m_3 ] with 
electron plasma frequency uo pe — y^Tm^, e _e 2 /m e = 5.64 x 10 7 s _1 x ^/ni2, e -[m -3 ], the cell size 
Ax = Ay = 0A4:c/(jj pe = 2.34 m/y/ n\2,e- [m -3 ] with 3x3 particles per cell and species, and the 
magnetic field amplitude |Bi2 1 = 8.94:m e c(jj pe /e = 2.9 mT x y / ni2, e -[ m_3 ]- The relativistic ion 
Larmor radius is defined as ru ~ c/uo c i — m p c 2 7i2/(e|Bi2|) = 1.9km x (m p / m e ) / yj ni2, e - [m -3 ] . 
Particles are symmetrically injected from both sides of the two-dimensional simulation box. 
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Fig. 4. — Simulation parameters 



3.1.1. Analysis of the particle spectra 

The comparison of the particle spectra is done at the same time in units of cj^ 1 = 6.3 x 
10 _6 s _1 x {m p / m e ) I yj n\2,e- [ni -3 ] . Figure [5] shows the spectra at t — 654 cj^ 1 . The electron 
spectra do not differ much from a thermal distribution, as well as the ion spectra, which have just 
thermalized, whereas the positron non-thermal tail is strong. A scaling of the maximum 7 with the 
ion to electron mass ratio is apparent. For the density ratio ni2,p+/^i2,e- — 0.7 the peak of the 
positron tail almost reaches the same level as the maximum of the thermal bulk and it is two orders 
of magnitude lower for ni2,p+/^i2,e- — 0-2. The comparison of the distributions of the different 
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species in Figure [5] (d) shows that the electron tail is much weaker than the positron tail. The ion 
spectrum has a completely different shape from the light species, which was also observed in the 
case of a pure electron-ion shock with high initial magnetization ( Sironi fc Spitkovsky||201l" ). 




Fig. 5. — Electron (a), positron (b) and ion (c) downstream spectra. Color codings are given in 
Fig. [4j (d) Comparison of electron (blue), positron (red) and ion (orange) spectra versus ~fm a /mi 
with m a /rrii the mass of the respective species normalized by the ion mass for a\ — 0.145 and 
mi / m e — 60. All spectra are plotted at t = 654 oo~^ . 

The spectra are fitted with functions of the form given in Equation Q, with the parameters 
listed in Fig. ?? in the Appendix, and used to calculate the jump conditions according to Equations 
([2]) and ([3]). In order to determine the adiabatic constant systematically, we also integrate the 
spectra numerically, with a deviation of the order of less than 0.1% deviation from the analytical 
result. The jump conditions are determined with the standard MHD model (S-MHD), where the 
adiabatic constant is 3/2, and compared to the advanced model (A-MHD) given by Equations ^ 
and The comparison with the simulation data is given in Fig. |6j By plotting the transversely 
averaged density against x\ and t the velocity of the moving shock front is measured, which is 
almost perfectly constant after a few 100's of oo~J~. To determine the density jump, the density 
was averaged over the full downstream region. This value is also constant during the entire shock 
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e volution. 

From Fig. [6] we see that the A-MHD model fits the simulation data better than the standard 
model, although the variations are only on the 1% level. The dependence of the shock speed on the 
magnetization o\ and ion kinetic energy U p /U e is in agreement with Figure |2j The variations for a 
constant total magnetization g\ and different energy ratios are rather small, but also here the trend 
towards higher density compression ratios and lower shock speeds for increasing ion kinetic energy 
ratio U p /U e is clearly recognizable. 
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Fig. 6. — Downstream parameters measured from simulation data and comparison of density com- 
pression and shock speed with theory (A-MHD), obtained from Equations and Q, and standard 
MHD theory (S-MHD) in brackets. 



3.1.2. Long time evolution of acceleration in mixed shock plasmas 

We also studied the temporal evolution of the adiabatic constant by integrating the spectra. 
The adiabatic constant of the positron spectrum initially decreases rapidly and is constant after 
a few oo~^ . The adiabatic constant of the electron spectrum increases first and then slowly drops 
towards the lower limit 3/2 of an ideal two-dimensional gas. The ion constant seems to follow the 
same trend, and at the end of our simulations, it is still in the increasing stage. 

Figure [8] (a)-(c) shows the temporal evolution of the maximum gamma of the non-thermal 
tail which finally determines the changes in the adiabatic constant. We find, that the maximum 
energy scales like 7 ma x oc (£ — to) a , with the values for a varying between ^ 1/3 and 1, as shown in 
Fig. [7| The positron energy increases faster than the electron energy, which is in agreement with 
the preferred energy transfer due to the synchrotron maser instability. We observe a scaling with 
the ion mass ratio rather than with the kinetic energy ratio. For a pure electron-ion shock in the 
unmagnetized case, a a = 0.6 has been observed (Fiuza et al. 2011, in preparation). The range 
of values for a obtained here are consistent with acceleration due to multiple scattering in small 



wavelength turbulence (determined by the collisionless length/time scales) as predicted by |Kirk fc 



Reville (2010). For Bohm diffusion, the spatial diffusion coefficient scales like k — \v/Z oc 



2 
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estimating the acceleration time ( Drury|[l983 ) 
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(8) 



as t acc oc 7 (Gargate & Spitkovsky 2011), where v Ul Vd are the upstream and downstream flow 
velocities and k u , the upstream and downstream spatial diffusion coefficients. According to |Kirk 



& Reville (2010), in the case of small-angle scattering the mean free path A is rather proportional 



to 7 , as well as the spatial diffusion coefficient, and therefore the maximum energy is expected to 
evolve as t 1 / 2 in the limit 7> 1. 

After t > 200 cj"- 1 the maximum positron energy stays almost constant (a < 0.1). 



Sironi & 



Spitkovsky (2011) analyzed the acceleration mechanisms in pure electron-ion plasmas in an oblique 
magnetic field with magnetization a\ = 0.1 and angle 9 — 75° to the longitudinal component for a 
mass ratio m p /m e = 16 and found the synchrotron maser instability to be the dominant process in 
such a configuration. The transverse electromagnetic wave modes affect mostly the electrons, which 
leads to a decrease of their longitudinal momentum, whereas the heavy ions propagate almost with 
the initial momentum. The electrons are accelerated towards the shock by the induced longitudinal 



wakefield (Lyubarsky 2006) and it was observed that both species enter the shock region with 



almost the same energy, so that the electric field does not persist in the downstream. In contrast, 
we observe a non-zero electric field in the downstream region (Fig. which is clearly above noise 
level. During the process of ion thermalization, with the characteristic spiral structure in the p2 —pi 



phase space (Hoshino & Arons 1991), the electric field is decreased due to the random motion of 



the ions, which causes the observed slowing down of the acceleration in Figure [8| The break in the 
power-law at t > 200 uo~^ appears when the ions have finally thermalized. The electric field in the 
far downstream region has reached its asymptotic value 0.25 | E12 1 at the same time. Moreover, at 
this stage, the positron energy becomes comparable to the ion energy, m e 72 5 e+/(^p72,p+) — 1, so 
that both positive species act in a similar way, whereas the electrons are still accelerated as their 
energy is small compared to the energy of the positive species. 
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Fig. 7. — a measured from the simulations according to 7 max oc (t — to) a . 



Furthermore, Sironi & Spitkovsky (2011) observed a decrease of the electron tail to a thermal 
spectrum for tu p i > 7000 for the above mentioned setup. We have not seen any indication of such a 
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200 400 600, 800 1000 1200 




Fig. 8. — Maximum energy evolution for electrons (a), positrons (b) and ions (c) with fits to a 
power-law of the type ^raax oc (t — to) a in (a) and (b) for electrons and positrons. The positron 
energies have been fitted with two power-laws. The index for low 7 is given in Fig. [7] together with 
the color codes. The index of the power-law for t > 200 uj~1 is less than 0.1. (d) shows the evolution 
of the downstream electric field normalized to the initial upstream field averaged over the entire 
downstream region (black), in the far downstream (red) and in the shock front region (blue) over a 
range of 100c/oj pe with exponentially decreasing fits. 

decrease as even the positron spectrum remains almost constant. The detailed analysis of the long- 
term evolution and the influence of the field structure on the acceleration rate for these scenarios 
will be discussed in a future work. 
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3.2. Decreasing the initial magnetic field 

We have performed a series of simulations with lower magnetizations than in the previous 
section in order to test the model for the jump conditions over a wide parameter range and compare 
them briefly to the case of a pure electron-positron plasma. The magnetic and electric fields are the 
same as in the previous setup. The two-dimensional simulation box consists of 5000 c/uj pe x 50 c/uj pe 
with cell size 0.2 c/u pe and 6 particles per cell and species. We also did tests with up to 25 particles 
per cell, showing negligible differences. The two beams interact at x\ — 300c/cj pe , which allows 
us to fully resolve the shock dynamics, but also to reduce the box size. The time step is again 
chosen as uo pe At = l/y/2 of the Courant condition in order to reduce simulation noise. In the run 
with m p /m e = 100 the interaction region is shifted to x\ — 1750 c/cu pei and the simulation box is 
increased to 10 4 c/uj pe in propagation direction and the total simulation time is 1.2 x lO 4 ^^ 1 . The 



different species configurations can be taken from Fig. 14 in the Appendix. 



3.2.1. Varying the ion fraction 

In the case of a pure pair plasma the electron and positron densities are nearly identical and 
show almost no spatial variation in the downstream. Only very weak filaments appear in the 
upstream region, which is an indication that the Weibel instability is almost completely suppressed. 
It is interesting to see that a magnetic precursor exists ahead of the shock also in the highly 
magnetized case. The phase spaces show a sharp transition between the downstream and the 
upstream region, also recognizable in the spatial density, which has been found to be a characteristic 
of superluminal shocks ( Sironi fc Spitkovsky||20lf ). The pair spectrum shows no evidence of a non- 



thermal tail and can be best fitted with / (7) = Cexp (— 7/A7) with A7 =17. 

Moving from pure pair plasmas to mixed configurations, we observe that the shock speed is 
decreased and the density compression ratio increased if the ion to electron mass ratio is increased, 
demonstrated in Figure [9] for a density ratio ni2, p +/ni2, e - = 0.2, which is in very good agree- 
ment with theory. The same behavior is found if the mass ratio is fixed and the density ratio 
is increased. The heavier ions cross the interaction region over a few electron skin depths, but 
soon they are reflected and the left-hand and right-hand populations are well-separated from each 
other. For demonstration purposes, the densities have been normalized to the upstream densities 
of each species, ^2,a/^i2,a? showing that all three particle components behave similarly. The two- 
dimensional density profiles reveal a weak filamentary structure for electrons and positrons, like in 
the case of the pair plasma, which is not apparent in the ion density. Similarly to the unmagnetized 
case in Stockem et al. (2011) the jump conditions only weakly depend on the real shape of the 
downstream spectra and the jump conditions in Fig. [14] are in good agreement with the theoret- 
ical model. The standard MHD model is a good approximation. Only at lower magnetizations 
we observe a stronger deviation from the simple model, in agreement with Figure [3j We observe 
that the final compression ratio is reached already after 3-4 kG 1 , showing a steady state after 20-30 
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Fig. 9. — Comparison of the averaged and normalized densities Ti2 a /nil a of electrons (blue), 
positrons (red) and ions (orange) with ni2,p+/ni2,e- — (a) and 0.2 else, m p /m e = 20 (b) and 
m p /m e = 100 (c). The simulation time is tuj pe = 3800. 



uo c f . |Haim et al. (2012) predict quasi-stationary solutions for low-a, but still not Weibel-governed, 
shocks only in case of electron-ion plasmas as the wave steepening will be stopped by energy dis- 
persion into whistler waves, which are not present in pair plasmas. Furthermore, if ions are present, 
the downstream structure along x\ shows strong wave generation on the scale of the ion Larmor 
radius. In the case of a pure electron-ion plasma, the structure is again much smoother and these 
oscillations become very weak, almost disappearing. The interaction between the three components 
is responsible for the oscillations present in the mixed component scenario, as previously observed 



in a ID setup by Hoshino et al. (1992). When the ions start to gyrate, the magnetic field is com- 



pressed, which leads to a compression of the pair density as well, because electrons and positrons 
are frozen into the field. The VB drift generates a current that reacts back to the magnetic field, 
reinforcing the downstream compressional oscillations. Indeed, we observe a perfect match between 
the out-of-plane magnetic field structure and the density of the pairs in the downstream. The ion 
density oscillates with the same frequency and a phase shift of 180°, reaching maximum density 
when the pair density is minimum. 



We observe that if the mass ratio is increased, the electron spectrum approaches the positron 
spectrum, which is due to the decreased total magnetization. But still in the run with a mass ratio 
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Fig. 10. — 2D (a) and ID (b) downstream distributions at tuj pe — 3800 for positrons (red), electrons 
(blue), ions (orange) with m p /m e = 20 and ni2,p+/^i2,e- = 0.6 and Maxwellian fit to the 2D data 
in black. 



of 100 and a x < 10~ 2 , the non-thermal electron tail stays weaker than the positron tail even for 
long simulation times. 

The analysis of the densities and magnetic fields showed that spatial variations along X2 are 
low. Also the differences in the distribution functions, obtained from ID and 2D simulations, are 



small (Figure 10), and arise essentially from the different statistics in ID vs. 2D simulations, which 
justifies to study the effects of a realistic mass ratio in ID simulations. 



3.2.2. Realistic mass ratio 

The ion mass ratio is further increased to a realistic proton to electron mass ratio m p /m e = 
1836, and the total magnetization is decreased to the limit a « 10 -3 , where the Weibel instability 



starts to become important (Spitkovsky 2005). Since the shock formation is determined by the 
proton cyclotron time scale uo~ p+ , we can study this process in detail. On the one hand, the reduced 
geometry was chosen due to limited computational resources, on the other hand, Weibel modes are 
suppressed and can be excluded as the shock driving mechanism. Although we are slightly above 
the threshold, we are aware that 2D effects might become important and plan to investigate their 
role in future work. 

Because of the large proton Larmor radius tl p+ = 8200c/cj pe for initial electron magnetization 
a e = 2, we increase the one-dimensional box to 4 x 10 4 c/cj pe with a cell size Ax± = 0.25c/cj pe , 
using 64 particles per cell and species. The interaction point of the counterpropagating beams is 
fixed at x\ — 10 4 c/(jj pe and the total simulation time is 4 x lO 4 ^"/ = 4.87 cj^Y with the proton 
cyclotron time scale uo~ p + ~ tl p +/c. Also in this case, the density ratios ni2,p+/^i2,e- — 0-2, 0.6, 1 
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are investigated. 

The runs with a realistic mass ratio show a highly dynamic structure in the beginning, which 
appears to be almost independent of the initial ion fraction. Figure 11_ shows the electron density, 
for the case ni2 } p+/ n i2 } e- — 0.6, against x\ and t which allows the determination of the shock 
velocity f3 = x\jt. Three stages have been identified from the analysis of the shock speed, which we 
discuss in detail for this density ratio. 




Fig. 11. — Electron density ri2^-/^i2,e- (&) and electron to proton density ratio 
[n2, e -/^i2,e-]/[^2,p+/^i2,p+] (b) for realistic mass ratio m p /m e = 1836 and ni2,p+/^i2,e- — 0-6. 
The dashed lines indicate x\jt — 0.49. 

In the first stage, which lasts until approximately 1.5 o;^, the ions are still cold and their 
phase space profiles differ much from that of the electrons and positrons. The counterpropagating 
beams are unaffected and propagate almost with the speed of light. At t « 4100 u~J- = 0.5 
the plasma is significantly compressed. At this stage we observe an extended ion gyro cycle, which 
reaches deeply into the downstream region, while the two populations of left and right electrons 
and positrons are almost separated and overlap only for a few electron gyro radii at the interaction 
region. The light particles are already thermalized and the beginning of a non-thermal profile is 
recognizable. While the electron distribution can be fitted well with a Maxwellian with thermal 
spread A7 =17, the positrons already show a clear non-thermal tail. 

In the second stage 1.5 < t < 3cj~p + , the ions start to respond slowly to the generated 
magnetic field compressed at the shock front and their distribution deviates from the initial cold 
upstream distribution, but the population is still far from being thermalized. The ion average 
density shows a strong spread around the electron and positron profiles and waves are generated 
on the scale of the proton Larmor radius. Electrons and positrons are further accelerated by the 



cyclotron instability (Amato & Arons 2006), with the strongest effect directly behind the shock 
front. In that region, the electron and positron spectra are almost equal, revealing a strong non- 
thermal component. The strong compression of the plasma particles directly behind the shock front 
is responsible for the rapid decrease of the shock velocity. A large amount of particles is reflected 
with the speed of light to both sides of the shock fronts at t « 12300 o^ 1 = 1.5 a;^. and the density 



structure in Figure 11 shows a second arc. The density profile becomes very dynamic, which has not 



been observed for low mass ratios and the phase spaces show large electron and positron momenta 
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where there is a strong mixture of the ion population, e. g. for tcj pe = 23800 at x\ — 1.4 x 10 4 c/uj pe . 

In the third stage, at approximately t « 24600 uj~} — 3 cj~ p+ , a quasi-steady structure is reached, 
where the shock speed matches the theoretical value /3 = 0.49. But even for these simulation sizes 
it is difficult to clearly identify the formation of the shock. The density profile in the downstream 
region oscillates in space on the ion scale c/cj c ^, revealing regions of proton or electron accumulation, 
as observed in the previous section for reduced mass ratios. These regions appear to be quasi-steady 



in time (Figure 11). 



The observed evolution of the shock speed is due to the involved dynamics of the different 
species. The electrostatic fields, which also appear in pure electron-ion shocks because of the 
different inertia, are only partially balanced, due to the presence of a light positive species. A 
precursor of electrons and positrons exists in front of the shock and the deceleration and acceleration 
of the light species resembles the crossing of the shock front in the Fermi acceleration process, which 
can enhance the acceleration of particles in mixed plasmas. 

In the other cases, ni2, p +/ni2, e - — 0.2 and 1, we observe a similar qualitative behavior with 
the appearance of three temporal stages. In all three cases, the electron spectra show only weak 
non-thermal acceleration and deviate most from a thermal spectrum for a low proton to electron 
density ratio (see Fig. fl2k), but the spread A7 of the peak energy of the thermal bulk increases 



with ni2 5j 9+/ni2,e-? as can be seen in Fig. 12i-c, and the highest tail energies are achieved in the 
case of a pure electron-proton plasma (see Fig. [12] :) . The maximum positron energy is independent 
of the proton fraction, but also here we observe an increase in the bulk spread. The proton spectra 
do not show evidence of non-thermal particle acceleration which is consistent with the results of 



Sironi & Spitkovsky (2011). 



3.3. Magnetic field in the plane 

We performed simulations with the upstream magnetic field in the plane for the total magneti- 
zation u\ — 0.145 and ion density ratio ni2,p+/ni2, e - = 0.7. We compare the cases 9 — 90°, 45°, 0° 
where 9 is the angle between the magnetic field and the longitudinal direction. The dominant ac- 
celeration process is determined by the magnetic field orientation, which was classified by |Sironi fc 



Spitkovsky (2011) into subluminal 9 < 9 cr i t ~ 34° and superluminal shocks 9 > 9 cr i t . Accordingly, 
particles gain energy in subluminal shocks by non-resonant interactions with Bell's waves and are 
efficiently accelerated while bouncing back and forth across the shock. In superluminal shocks, if 



ions are present, the synchrotron maser waves transfer energy to the electrons. jSironi fc Spitkovsky 



(2011) observed a short power-law tail stemming rather from heating than acceleration. 



Our results are in agreement with these findings for the electrons. At the end of the simulation, 
at tuj pe = 6000, the electron spectra show no sign of non-thermal acceleration if the magnetic field 



is superluminal, see Figure 13 In the subluminal case, the highest energies are achieved, which 



reach the level of that of the positrons. However, the positrons show a different behaviour. The 
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Fig. 12. — Electron (blue), positron (red) and proton (orange) distributions versus ^ym a /m p with 
m a /m p the mass of the respective species normalized by the proton mass. The spectra have been 
averaged over the entire downstream region at tuj pe — 4 x 10 4 . The initial density ratios are 
ni2,p+/ni2,e- = 0.2 (a), 0.6 (b) and 1 (c). 



maximum positron energy is the same in all three runs, with the strongest tail (=higher fraction of 
energy in the non-thermal particles) in the case of the perpendicular shock. 

The ion spectra reach the same maximum kinetic energy as the positrons. The shape of the 
particle distributions varies much more than for the other components and resembles a thermal 
distribution only in the case of a perpendicular initial magnetic field. For the cases with a par- 
allel component, the spectra become narrow, rather like a 3D-Maxwellian as it was observed for 



superluminal shocks also in Sironi & Spitkovsky (2011). 




Y 

Fig. 13. — (a) Electron, (b) positron and (c) ion distributions for tuj pe — 6000 and angle between 
upstream magnetic field and longitudinal direction 9 = 90° (blue), 9 — 45° (green), = 0° (brown). 

4. Discussion 

We have investigated the shock generation in plasmas consisting of electrons, positrons and 
protons for different perpendicular upstream magnetic fields. The standard one-fluid jump con- 
ditions have been extended for a multi-component plasma and the real shape of the downstream 
particle distributions has been taken into account. The calculations predict a decrease of the shock 
speed if either the mass ratio or the upstream ion fraction are increased. The higher the upstream 
magnetization, the less important become the effects of the real particle distribution. If the devi- 
ations from a Maxwellian are low, or if the magnetization is high enough, for a highly relativistic 
upstream the shock speed is determined by a simple second order equation, which depends only on 
the effective upstream magnetization. 

Simulations of shocks in mixed plasmas have been performed for a constant total magnetiza- 



- 19 - 



tion with different initial ion kinetic energies and compared with our advanced theoretical model. 
The shock has been launched by injecting two counterpropagating beams from each side of the 
two-dimensional simulation box. Whereas the standard model predicts the same jump conditions 
independent of the ion kinetic energy, the advanced model fits better and predicts the increase in the 
compression ratio and decrease of the shock speed for increasing ion kinetic energy. Nevertheless, 
the differences are on the 1% level only. 

The evolution of the maximum energy was found to be a power-law with a power 1/3 < a < 1, 
which depends on the initial parameters, but clearly indicates that the acceleration process is slower 
than the usually considered Bohm diffusion and is consistent with scattering off small-scale magnetic 
fluctuations. The non-thermal tail of the positrons was found to be constant after tuo~^ « 200, which 
coincides with the thermalization of the ions. 

The theoretical model has been tested for a wide range of parameters, showing a good agreement 
between simulation and theory, but a weak dependence on the actual shape of the spectra. If the 
ion fraction is increased, the evolution of the density is similar in all components on time scales of 
the inverse ion cyclotron frequency. During this transition phase the light particles run ahead of the 
heavy ion species, undergoing an oscillation along the shock propagation direction. This separation 
of the species has been observed for a mixed plasma even in the case of no initial magnetization, 
whereas in a pure electron-ion plasma both species are always perfectly matched right from the 
beginning of the shock formation like in the case of a pure pair plasma. 

We observed that two-dimensional effects become more important if the total magnetization is 
low due to the increased Larmor radius. Nevertheless, transverse spatial dependences were found 
to be low and the particle spectra are similar. 

This work was partially supported by the European Research Council (ERC-2010-AdG Grant 
267841) and FCT (Portugal) grants SFRH/BPD/65008/2009, SFRH/BD/38952/2007, and PTDC/FIS/111720/200 
Simulations were performed at the 1ST cluster (Lisbon, Portugal). 



A. Derivation of the jump conditions 



The purpose of this section is to derive an expression for the shock speed in terms of known 
upstream quantities, with which the jump conditions for a scenario described in Figure [T] can be 
determined. The following calculations base on the conservation equations for a single fluid in the 



paper by Kennel & Coroniti (1984), which are expressed in the shock frame, where upstream and 



downstream components both propagate perpendicular to the shock front. The magnetic field is 
oriented perpendicular to the shock front, as well as to the propagation direction of the particles, 
as demonstrated in Figure [Tj The conservation equations in the shock frame for multiple species 
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are thus given by 



n2,aU2s 



PlsBls — fastis 

a 2 i Pl-sBi s 2 i #2s-B|s 
PUlu 2^ n haHl,a + — T = P2sl2s n ^2,a + — T 



47T 

B 



(Al) 
(A2) 

(A3) 
(A4) 



with proper velocity U{ s = Pisjis, perpendicular magnetic field component Bi Sl specific enthalpy 
Mi, a = ^ a c 2 + (ei ia + Pi, a ) /^i,a: rest frame energy density e^ a and pressure density p^ a . In the 
following;, we assume a cold upstream with p\ a — c\ a = 0. Similar to the case of a single fluid, a 
downstream adiabatic constant T a can be defined for each species by the pressure-energy relation 
Pi,a = (T a - l)(ei, a - n^ a m a c 2 ). Thus, the sum 
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_ p 6 — — . Combining Equations (|Al|)-(|A4|) yields the determination equation of the shock 
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speed in the shock frame, given by 
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with Y :— Pu/p\s an d w i '■— 1 + ^ 2 ' e+ + ^ 2 ' p+ ■ The total magnetization is defined as 
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(A7) 



according to Hoshino et al. (1992), providing the total upstream magnetization a± = B\ 2 j (47mi2, e -Wi77i e c 2 7i2). 
Performing a Lorentz transformation into the downstream frame, using n^ a — Tiij^/jij^ &i — 
Bijhij, := 0*2 = -fcs > 0, 7 := 72, = 7,2, Pu = -(\PM + P)/( l + \Pu\P) and lu = 
7712(1 + 0|0i2|) one obtains the determination equation for the shock speed with parameters de- 
fined in the downstream ( = simulation) frame 
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which is an algebraic equation of fifth order in /3. It is only slightly simplified by the assumption of 
equal downstream densities, applying W2/W1 = 1. For a highly relativistic approximation 712 ^> 1, 



Equation (A8) reduces to a quadratic equation in the shock speed 
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Fig. 14. — Downstream parameters measured from simulation data and comparison of density 
compression and shock speed with theory in brackets, obtained from Equations (fTl) and ([3]). 
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